Sphere penetration by impact in a granular medium: 
A collisional process 
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Abstract. - The penetration by a gravity driven impact of a solid sphere into a granular medium 
is studied by two-dimensional simulations. The scaling laws observed experimentally for both 
the final penetration depth and the stopping time with the relevant physical parameters are 
here recovered numerically without the consideration of any microscopic solid friction but with 
dissipative collisions only. Dissipative collisional processes are thus found as essential in catching 
the penetration dynamics in granular matter whereas microscopic frictional processes can only be 
considered as secondary effects. 
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Introduction. — A better knowledge of the impact of 
a solid object into a granular target has many applied in- 
terests from both geologist and ballistic point of view [1,2]. 
Since a few years, numerous studies have been carried 
out by physicists interested in the ejection process [3-8], 
the crater morphology [9-12] and the penetration dynam- 
ics [4, 13-20], searching for scaling laws for the crater size 
[9-12] and penetration depth [13-23]. For the ejection pro- 
cess, a spectacular thin granular jet raising very high can 
be observed after the impact on small grains of rather low 
packing fraction [3-7] whereas an opening granular corona 
is seen for larger grains of rather high packing fraction [8] : 
the effect of air is crucial in the former case [5,6] whereas it 
is negligible in the latter one. For the crater morphology, 
the scaling laws found for high energy impacts of plan- 
etary interest [1] stand for low energy impacts of small 
scale laboratory experiments [9-12], indicating some uni- 
versal physical processes involved in the crater formation. 
For the penetration dynamics, the observed deceleration 
of the impacting sphere towards its final stop is usually ex- 
plained by a complex drag force resulting from frictional 
and collisional processes and involving several terms: a 
linear depth dependent term [17] arising from solid fric- 
tion and velocity-dependent terms of linear or quadratic 
form arising from the ballistic [14, 17, 18] or the fluid me- 
chanics point of view [21]. Such different force terms 
are not easy to extract from the sphere trajectory usu- 



ally tracked by video means [13-17], but recently direct 
force measurements with an accelerometer inside the im- 
pacting sphere [19] reveal that a force proportional to the 
velocity squared is indeed experienced by the impacting 
sphere at least during its first penetration stage at high 
velocity and shallow depth. The stopping time of the 
sphere has been studied in different experimental works 
and displays a striking and rather counter-intuitive be- 
haviour: it is a decreasing function of the impact velocity 
with an asymptotic plateau for large enough impact ve- 
locities [13,17,19]. The characteristic time scaling for the 
plateau value was proposed to be r = (d/g) 1 ^ 2 [13,17] or 
r = {p/ p g ) 1 / i {d/ g) 1 / 2 [19] for a velocity larger than the 
typical characteristic velocity V = {gd) 1 / 2 [17,19], where 
p g is the density of the grains, d and p are the diameter 
and the density of the sphere. The observed scaling law 
for the final penetration depth S that may be written as 
S/d oc (p/p g )P(H/d) a , where H is the total falling dis- 
tance covered from release to rest, is not yet satisfactorily 
explained, as well as the values of the two power expo- 
nents a and /3 found to be around 1/2 in experimental or 
numerical works [10,14,18,20-22]. Experiments are es- 
sentially 3D (except the real 2D experiment of Ciamarra 
et al. [13]) whereas the numerical simulations are 2D with 
similar values for the power exponents. The existence of 
a finite penetration depth S is always referred to the ex- 
istence of a non-zero solid friction p in contrast to the 
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Fig. 1: Time evolution of the penetration z of a sphere im- 
pacting the granular medium, where z = corresponds to the 
initial surface of the grain piling. The projectile of diameter 
d = 30 d g = 30 mm and of density p — p g — 2520 kg m~ 3 hits 
the grains at the velocity v t = 1ms -1 . Inset: snapshot during 
penetration by impact. 

case of simple fluids where the sphere would not stop but 
reach a limit velocity in the absence of solid friction. The 
(j, dependence of S was proposed to be 6 oc fi^ 1 [10] (the 
penetration depth would thus be no more finite for zero 
solid friction) from experimental investigation varying the 
grain material. But experimentally, it is hard to change 
the coefficient of solid friction in a large range and without 
changing other crucial parameters such as the coefficient 
of restitution and the solid fraction. 

In the present paper, we show by numerical simula- 
tions that no microscopic solid friction is necessary to 
explain a finite penetration depth for a sphere impact- 
ing on granular matter. Furthermore, we show that the 
scaling law 8/d oc [pj pg) 13 (H/d) a observed experimen- 
tally or numerically with non-zero microscopic solid fric- 
tion still stands in the zero microscopic solid friction limit. 
We also recover that the stopping time t s is constant at 
large enough impact velocities and show that it scales as 
t s ~ t = (pd/ pgg) 1 / 2 for impact velocities larger that the 
typical characteristic velocity V = [pgdj pg) 1 ! 2 . These 
scalings clarify the previous scalings discussed above. The 
numerical results will be shown and discussed after hav- 
ing presented the 2D numerical method used in the present 
paper. 

Numerical method. — We use the method of molec- 
ular dynamics to perform two-dimensional simulations in 
the geometry shown in the inset of fig. Q] The granular 
target is prepared by the sedimentation of an initial dilute 
configuration under the action of the gravity acceleration 
g = 9.8 ms . The grains are modeled as a random pack- 
ing of spheres of mean diameter d g , mass m g and density 
p g contained in a rectangular box bounded by hard walls. 
In order to avoid any crystallization of the packing, a slight 
dispersion in the grain diameter is introduced, with an uni- 



form distribution in the range 0.8d g to 1.2d g . Before the 
impact, the packing fraction in the granular medium is 
<f> ~ 0.83. The projectile is a sphere of diameter d, mass 
m and density p which is dropped onto the granular pack- 
ing. The projectile is thrown downwards at the velocity 
required, and its value at impact will be noted Uj. This 
is equivalent to the usual experimental situation where 
the sphere is dropped from the height h and impacts the 
granular surface with the velocity V{ = \J2gh. Note that 
the box containing the granular medium is large enough 
(> 8d) so that the projectile is not affected by the sur- 
rounding walls and the layer of grains is high enough to 
avoid any bottom wall effects during the penetration [22]. 
The number of grains in the simulations ranges from 10 
for small boxes to 10 5 for largest ones. 

As the goal of this paper is to show that microscopic 
solid friction is not essential in explaining the finite pene- 
tration of a projectile into a granular material, we do not 
take into account any static nor dynamic friction between 
the grains. The interaction forces are thus taken as purely 
normal with no tangential components. The interaction 
force between two grains, or between one grain and the 
projectile or the bounded walls, is modeled as a dissipa- 
tivc Hertz law [24] such as 

F n = - 7 §, (1) 

where k is the non-linear stiffness, 7 is a damping coeffi- 
cient, and £ and d^/dt are respectively the interpenetra- 
tion and the velocity of interpenetration of the grains. For 
two identical spherical grains of diameter d g , Young mod- 
ulus E and Poisson coefficient v, the non-linear stiffness k 
is given by k = 2Ey/2d g /3(l — v 2 ). The collision time for 
a non-dissipative contact (7 = 0) between two grains is 

r c « 3.21 (2=2) v ~y^ ( 2 ) 

where v n is the relative normal velocity and m ff = m g /2 
the effective mass for two identical colliding grains. The 
collision time is not very different in the case of a non- 
zero dissipation (7 ^ 0). The equations of motion for 
the grains are integrated using a standard second order 
Verlet algorithm. The choice for the numerical time step 
At must be such that At -C t c in order to ensure numerical 
accuracy. 

In the following, we present numerical simulations for 
a granular material composed of glass spheres (density 
Pg = 2520 kg m -3 ) with an effective elastic modulus E* = 
E/(l — v 2 ) = 69 x 10 9 Pa and of mean diameter d g = 1 mm 
(mass m g = 1.3 mg). The non-linear stiffness is thus 
k = 2 x 10 9 Nm 3 / 2 and the collision time for a typical 
collision velocity v n = 1ms -1 is t c = 2.7 /is. The time 
step is chosen as At = 0.1 fis <C t c . For a non-zero 
damping coefficient 7, the coefficient of normal restitu- 
tion for normal incidence e n = —v^/v^, which is the ratio 
of the relative normal velocity after an impact over the 
velocity before the impact, is smaller than the ideal limit 
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Fig. 2: Normalized penetration depth S/d of the sphere as 
a function of the normalized total falling distance H/d for 
different projectile diameters: d = 20d 9 (V), d = 30d 9 (o), 
d = 40d 9 (□) and d = 70d 9 (A). The sphere/grains density 

ratio is p/p g = 1. ( ) Power law fit S/d oc (H/d) a , with 

a ~ 0.31. 



value 1 for perfect elastic collisions. More precisely, with 
model cq. (JTJ , e„ decreases slowly with the collision veloc- 
ity as (1 — e„) oc (u*) -1 / 5 [24]. With the damping value 
7 = 0.065kgs _1 , the restitution coefficient is e„ = 0.9 
for v l n = lms -1 . It should be stressed that the follow- 
ing results are qualitatively only weakly dependent on the 
damping factor value. 

Numerical results and discussion. — Figure [T] dis- 
plays the time evolution of the position z(t) of the pro- 
jectile during its penetration through the grains, where z 
is the distance between the initial horizontal free surface 
of the granular medium and the bottom of the impacting 
sphere. The penetration increases with time up to a max- 
imum depth S at the time t s . Note that a small rise of 
the sphere, of typically a few percent of the total pene- 
tration, is observed at the end of most of the runs. We 
attribute this effect to the absence of microscopic static 
friction in the interaction law [eq. (pQ)]. Indeed, the grains 
ejected during the collision process are redeposited on the 
granular material and exert an increasing pressure on it, 
inducing a downwards motion of the granular material. 
This small effect was not reported in numerical simula- 
tions including microscopic solid friction [18]. In the fol- 
lowing, we shall consider this maximum penetration 5 as 
the final penetration depth and the corresponding time t s 
as the stopping time. Note that dropping the sphere from 
slightly different x-positions gives very similar z(t) even if 
the acceleration signals are quite different. 

Penetration depth. Figure [2] shows the evolution of 
the normalized penetration depth S/d as a function of the 
normalized total falling distance H/d, where H = h + 
S is the sum of the free-fall height h = vf /2g and the 
penetration depth. The diameter of the impacting sphere 
ranges from d = 20d g to 70d g , and its density is here 



Fig. 3: (5/d)(H/d) a as a function of the sphere/grains density 
ratio p/p g . (- -) Power law fit (8/d)(H/d)- a oc (p/p g ) p , with 
(3 ~ 0.74. 

kept constant and equalled to the grain density (p = p g ). 
In the log-log plot of fig. the data are aligned along a 
straight line of slope a = 0.31 ± 0.02 meaning that the 
penetration depth 8 varies with H following a power law 
6/d oc (H/d) a . The variation of the penetration depth 
with the density (i.e. with the mass) of the projectile 
is shown in fig. [3J where (5 / d) / (H / d)~ a is plotted as a 
function of the density ratio p/p g between the projectile 
and the grains. The penetration depth is found to depend 
on p/ p g with a power law of the form S/d oc (p/pg) 13 , with 
[3 = 0.74 ± 0.02. Finally, grouping the fall height and the 
mass dependencies, one obtains 



d \ p„ ) V d 



(3) 



with A = 0.92 ± 0.02. This power law scaling for the 
penetration depth is observed in various experimental and 
numerical studies [10, 14, 18, 20-22]. The value of the ex- 
ponent a ~ 0.31 is close to the commonly reported values 
between 0.33 and 0.40, and the exponent f3 ~ 0.74 char- 
acterizing the dependence with the density ratio is not 
far from the reported values ranging from 0.50 to 0.61. 
Thus, the impacting sphere stops at a finite depth with- 
out any microscopic friction and the scaling laws for the 
penetration depth are in agreement with those observed 
experimentally or numerically with friction. 

Forces. Let us now examine the forces undergone by 
the sphere after the impact, to extract the main ingredi- 
ents leading to its stop. The force exerted by the grains 
on the sphere depends mainly on the penetration depth z 
and on the velocity v of the projectile [17-19]. For dense 
packings, experimental results show that this resistance 
force may be separated into two independent functions of 
position F z (z) and velocity F v (v) so that the Newton law 
for the sphere motion can be written as 



dv 
"dT 



= fj 



F z (z) F v (v) 



(4) 
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Fig. 4: Reduced acceleration — du/df + g of the impacting 
sphere (diameter d — 20d g , density p — p g ) as a function of 
its velocity squared v 2 , at a specific depth: (o) z = 30 mm, 
(□) z — 40 mm and (A) z = 50 mm. ( — ) Linear trends of the 
data by — dv/dt + g — F z {z)/m + v 2 /d\, with d\ ~ 36.4mm. 



F z is usually attributed to solid friction [4, 17] and F v is 
considered of collisional or inertial origin [17-19]. 

Figure 2] displays the reduced acceleration of the sphere 
as a function of its velocity at a given penetration depth 
(z = 30 mm) and suggests that the drag force F v is pro- 
portional to v 2 , in agreement with the expression mv 2 jd\ 
proposed recently by Katsuragi and Durian [17], where d\ 
is a characteristic dissipative length. Note that the re- 
duced force is non-zero at vanishing velocity as there is 
still a non-zero depth dependent force term F z (z). This 
behavior is observed with different sphere diameters in the 
range 20 mm < d < 100 mm and different density ratio in 
the range 1 < p/p g < 10. Figured shows that d\ does not 
depend on the depth of the projectile z and thus F v does 
not depend on z as found experimentally [17], which jus- 
tify the writing of eq. ((H). Figures [5p and[SJ; show that d\ 
is proportional both to the projectile diameter d and to the 
density ratio pj p g leading to d\ cx pdj p g . The velocity de- 
pendent force term F v (v) thus scales as F v /m cx p g v 2 / pd 
indicating its inertial or collisional origin. 

Plotting now in fig.[fjh the force term at vanishing veloc- 
ity uasa function of z shows that the simple dependence 
F z cx z proposed by Katsuragi and Durian [17] is com- 
patible with our data despite the high scattering at low 
z i z ^ d) where the sphere is not fully immersed in the 
granular medium. Figures and(BJi show that F z /m is 
proportional to 1/d and p g /p. The depth dependent force 
term F z thus scales as F z (z)/m cx p g gz/pd. This force 
term linear in z has been previously proposed and seen by 
various authors [17-19] with a solid friction origin, but the 
extracted coefficient of friction necessary to fit the data 
was far from the standard values [17]. Here, the F z {z) 
force term does not come from microscopic solid friction 
as there is no microscopic solid friction in our numerical 
simulations. The depth dependent term F z (z) can thus 




Fig. 5: (a) Characteristic length di as a function of the depth 
z of the projectile. ( — ) Constant fit di = 36.4 mm. (b) Char- 
acteristic length di as a function of the diameter d of the pro- 
jectile (density p = p g ) at a depth z = 20 mm. ( — ) Linear fit 
di = 1.3d. (c) Characteristic length di/d as a function of the 
density ratio p/pg at z = 30 mm for a projectile of diameter 
d = 20d g . ( ) Linear fit di/d = 1.2p/p g . 



simply viewed as a hydrostatic term as already suggested 
in [25] for the experimental penetration of flat plates and 
in [4] for an impacting sphere in very loose sand. The pres- 
sure increases linearly with the depth as p g gz and so is the 
force F z on the sphere. This "hydrostatic" force term is 
not Archimedean as the penetrating sphere is never im- 
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the projectile during its penetration may be written as 



Fig. 6: (a) F z /m as a function of the depth z of the projectile 
(diameter d — 20d g , density p — p g ); (b) F z /m as a function 
of the inverse of the diameter 1/d of the projectile (density 
p = p g ) at a depth z = 20 mm; (c) F z /mg as a function of 
the inverse of the density ratio (p/p 9 ) _1 at z = 30 mm for a 

projectile of diameter d = 2Qd g . ( ) Linear fits of numerical 

data. 



merscd in the granular packing before it stops as can be 
seen here numerically (see the inset snapshot of fig.Q} and 
experimentally [19]. It follows that the Newton's law for 



dt> 



Pg9 

O-z — T z 

pa 



PgV_ 

pd 



(5) 



where a z = 1 ± 0.1 and a v = 0.8 ± 0.1 are numerical 
prefactors. 

Stopping time. Let us now focus on the stopping time. 
The inset of fig. [7] displays the stopping time t s of the pro- 
jectile as a function of its impact velocity Vi. Numerical 
data are observed to be very scattered depending on the 
velocity, diameter and density of the projectile. More pre- 
cisely, t s decreases with the impact velocity Vi for any 
given size and density of the projectile and increases with 
the projectile diameter d and with its density ratio p/p g - 
Two time scales can be extracted from the non-linear dif- 
ferential equation ([5|) by considering independently the 
two force terms. Considering only the depth dependent 
force term, the characteristic time t z for the stopping time 

,Pg9J 

which is independent of the impact velocity. Considering 
now only the non-linear velocity dependent force term in 
cq. (|5|) leads to the characteristic time 



pd 

PgVt 



(7) 



which depends on the impact velocity Vi. As both force 
terms play a non-negligible role in the penetration (the 
velocity dependent force term decreases from its max- 
imal value at impact to zero at the stop whereas the 
depth dependent force term increases from zero at im- 
pact to its maximal value at the stop), the total stop- 
ping time t s (T z ,T v ) is a combination of the two charac- 
teristic time scales t z and t v , and can thus be expressed 
as t s jr z = f(r z /T v ). Note that the time scale ratio is 
t z /t v — Vi(p g /pgd) 1 ^ 2 , which corresponds also to the ve- 
locity ratio Vi/V of the impact velocity Vi to the charac- 
teristic velocity V = (pgd/ Pg) 1 ^ 2 ■ 

By using the characteristic time scale r z and velocity 
scale V, i.e. by plotting t s /r z as a function of Vi/V, all 
the data collapse on the same master curve with two dis- 
tinct parts (fig. [7|): (i) for low enough impact velocities 
{ v ilV ^ 2), the stopping time decreases with increasing 
impact velocity; (ii) for large enough impact velocities 
( V i/V 2), a plateau is reached, so that the stopping 
time docs not depend on the impact velocity and tends to 
the constant value t s /r z ~ 1.7. The existence of a single 
curve in the normalized plot t s /r z versus Vi/V, and the 
critical values Vi/V ~ 2 and t s /r z ~ 1.7 close to 1 indicate 
that the typical velocity scale V and time scale t z are rel- 
evant physical parameters for the problem, which validate 
the two independent model forces acting on the penetrat- 
ing sphere. The stopping time calculated by eq. (0) with 
a stop condition at v = agrees with the simulation data 
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Fig. 7: Normalized stopping time t a jr as a function of the 
normalized impact velocity Vi/V, for: p/p g = 0.5 and (▼) d = 
20mm, (■) d = 30mm, (■*) d = 40mm, (►) d = 60mm, 
(•) d = 80 mm, ( A) d = 100 mm; /9//9 s = 1 and ( V) d = 20 mm, 
(□) d = 30 mm, (<a) d = 40 mm, (o) d = 70 mm; p/p g = 1.5 and 
(♦) d = 30 mm; p/pg = 2 and (>) d = 20 mm, (A) d = 30 mm; 
p/pg — 3 and (o) d = 20 mm; p/p s = 4 and (+) d = 30 mm; 
p/p g = 6 and (x) d = 20 mm. ( — ) Guideline for the eyes. ( — 
) Calculated stopping time from model eq. ([5}. Inset: stopping 
time t s as a function of the impact velocity Vi for the same set 
of data. 



to within 20% (fig. [7]). The calculated penetration depth 
S/d obtained by solving eq. ([5]) deviates from the data by 
as much as 50% since small errors in the force terms in the 
approximated model equation have a more important ef- 
fect on the integrated position than on the stopping time. 

Conclusion. — By performing numerical simulations 
for an impacting sphere in a frictionlcss granular mate- 
rial, one obtains the same scaling law for the penetration 
depth as in simulations with solid friction or real exper- 
iments. This shows that dissipation by microscopic solid 
friction can be ignored and dissipation by collisions is suf- 
ficient to catch the penetration dynamics in granular mat- 
ter. Analysing the forces reveals that besides a velocity 
squared force term from collisional origin, exists a depth 
dependent force term. This pressure increasing term with 
the depth is sufficient to stop the sphere, and steric hin- 
drance with dissipation prevents the sphere from settling, 
by contrast to simple liquids. The scalings found for the 
two force terms allow for the prediction of the scaling of 
the stopping time which is indeed observed: a plateau 
value at high impact velocity and an increasing value at 
decreasing velocity. This scaling for the stopping time 
without any microscopic solid friction is also in agreement 
with what was observed previously. If microscopic solid 
friction appears needless to explain the scaling laws for 
both the final penetration depth and the stopping time of 
the projectile, it would certainly affect their precise value 
from a quantitative point of view. 
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